Read in data and results

gene_count_D9_df_filtered <- readRDS(here("bulkRNAseq", "results", "gene_count_D9_df_filtered.rds"))
gene_count_D30_df_filtered <- readRDS(here("bulkRNAseq", "results", "gene_count_D30_df_filtered.rds"))
sample_D9_metadata_df <- read_tsv(here("bulkRNAseq", "results", "sample_metadata_D9_df.txt"))
## Rows: 23 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): sample, day, mice, genotype, exh_subset, sample_num
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
sample_D30_metadata_df <- read_tsv(here("bulkRNAseq", "results", "sample_metadata_D30_df.txt"))
## Rows: 12 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): sample, day, mice, genotype, exh_subset, sample_num
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
norm_counts_D9 <- readRDS(here("bulkRNAseq", "results", "deseq_norm_counts_D9_exh.prog.vs.term.rds")) %>% 
  as_tibble(rownames = "gene")

norm_counts_D30 <- readRDS(here("bulkRNAseq", "results", "deseq_norm_counts_D30_exh.prog.vs.term.rds")) %>% 
  as_tibble(rownames = "gene")
dge_D30_compare_exh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_exh.prog.vs.term.txt")) # WT
## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D9_compare_exh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_exh.prog.vs.term.txt")) # WT
## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D9_progexh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_progexh_genotype.Enhdel.v.WT.txt"))
## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D9_transexh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_transexh_genotype.Enhdel.v.WT.txt"))
## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D9_termexh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_termexh_genotype.Enhdel.v.WT.txt"))
## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D30_progexh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_progexh_genotype.Enhdel.v.WT.txt"))
## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D30_transexh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_transexh_genotype.Enhdel.v.WT.txt"))
## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dge_D30_termexh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_termexh_genotype.Enhdel.v.WT.txt"))
## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
sc_dge <- read_tsv(here("scRNAseq", "data", "processed_data_tables", "01_dge.genotype_one-v-one_per-cluster.tsv"))
## Rows: 216930 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): gene, contrast, Clusters
## dbl (5): p_val, avg_logFC, pct.1, pct.2, p_val_adj
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gs_df <- read_tsv(here("scRNAseq", "data", "gene_signatures",
                       "02_gene_signatures.msigdbr_H-C2-C7.tsv"))
## Rows: 1509466 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (13): gs_cat, gs_subcat, gs_name, gene_symbol, ensembl_gene, human_gene_...
## dbl  (5): entrez_gene, human_entrez_gene, gs_pmid, taxon_id, num_ortholog_so...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gs_subset_bp_df <- gs_df %>%
  filter(str_detect(gs_name, "HALLMARK|REACTOME"))
length(unique(gs_subset_bp_df$gs_name))
## [1] 1665
# 1665
gs_subset_bp_list <- gs_subset_bp_df %>%
  named_group_split(gs_name) %>%
  map(~..1$gene_symbol)

gs_subset_tcell_df <- gs_df %>%
  filter(str_detect(gs_name, "_T_CELL|TCELL|CD8|CD4|EXH")) %>%
  filter(str_detect(gs_name, "HALLMARK|REACTOME", negate = T)) %>%
  filter(str_detect(gs_name, "NEUTROPHIL|MONO|ERYTHROBLAST|DC|NKT|LIN|BCELL|LSK|CD40|CD44|MAST", negate = T))
length(unique(gs_subset_tcell_df$gs_name))
## [1] 1446
# 1446
gs_subset_tcell_list <- gs_subset_tcell_df %>%
  named_group_split(gs_name) %>%
  map(~..1$gene_symbol)
gsea_D9_tcell <- read_tsv(
  here("bulkRNAseq", "results", "gsea_tcell_wide_D9_genotype.Enhdel.v.WT.txt")
)
## Rows: 1445 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr  (1): gs_name
## dbl (12): p.val_progexh, p.val_termexh, p.val_transexh, q.val_progexh, q.val...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gsea_D9_bp <- read_tsv(
  here("bulkRNAseq", "results", "gsea_bp_wide_D9_genotype.Enhdel.v.WT.txt")
)
## Rows: 1660 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr  (1): gs_name
## dbl (12): p.val_progexh, p.val_termexh, p.val_transexh, q.val_progexh, q.val...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

DGE, term vs prog exh

my_plot_volcano <- function(dge_df, label_genes) {
  dge_df %>% 
    drop_na(padj) %>% 
    ggplot() + 
    aes(log2FoldChange, -log10(padj)) + 
    geom_point(size = 0.5, alpha = 0.5) + 
    geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "red") + 
    geom_point(aes(color = padj < 0.05), 
               data = ~..1 %>% filter(gene %in% c(label_genes)), 
               size = 2) + 
    geom_text_repel(aes(label = gene, color = padj < 0.05), 
                    data = ~..1 %>% filter(gene %in% c(label_genes)), 
                    max.overlaps = Inf) + 
    scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) + 
    theme(legend.position = "none")
}
label_genes <- c("Havcr2", "Slamf6")
dge_D30_compare_exh %>% 
  my_plot_volcano(label_genes = label_genes) + 
  labs(title= "D30")

dge_D9_compare_exh %>% 
  my_plot_volcano(label_genes = label_genes) + 
  labs(title= "D9")

# gene list from miller et al 2019
inh_receptors <- c("Pdcd1", "Havcr2", "Lag3", "Cd244", "Entpd1", "Cd38", "Cd101", "Tigit", "Ctla4")
surface_receptors <- c("Tnfsf9", "Tnfrsf4", "Il2rb", "Klrg1", "Cd28", "Icos", "Tnfsf14", "Sell", "Il7r")
eff_molecules <- c("Ifng", "Il10", "Gzma", "Gzmb", "Prf1", "Tnfsf10", "Fasl", "Tnf", "Il2")
tfs <- c("Nfkb2", "Id2", "Hif1a", "Bach2", "Nfkb1", "Id3", "Tcf7", "Lef1", "Satb1", "Batf", 
         "Nfatc1", "Prdm1", "Runx1", "Runx3", "Tbx21", "Tox", "Bcl6", "Eomes")
chemokines <- c("Cx3cr1", "Ccl5", "Ccl4", "Ccl3", "Csf1", "Cxcr5", "Ccr7", "Xcl1", "Cxcl10")
select_genes <- c(inh_receptors, surface_receptors, eff_molecules, tfs, chemokines)
plot_gene_heatmap <- function(norm_counts, gene_list, 
                              col_index_ordering = c(3, 6, 9, 11, 2, 5, 8, 10, 1, 4, 7), 
                              hheight = unit(1.5, "in"), 
                              hwidth = unit(2, "in")) {
  h_mtx <-
    norm_counts %>% 
    filter(gene %in% gene_list) %>% 
    column_to_rownames("gene") %>% 
    as.matrix() %>% 
    t() %>% 
    scale() %>% 
    t() %>% 
    .[rownames(.)[order(match(rownames(.), gene_list))], 
      col_index_ordering]
  Heatmap(
    h_mtx,
    width = hwidth, 
    height = hheight, 
    col = colorRamp2(colors = c("purple4","yellow"), 
                       breaks = c(-3, 3)), 
    name = "row scaled\nnorm counts", 
    cluster_rows = FALSE, 
    cluster_columns = FALSE
  )
}
plot_gene_heatmap(norm_counts_D9, inh_receptors)

plot_gene_heatmap(norm_counts_D9, surface_receptors)

plot_gene_heatmap(norm_counts_D9, eff_molecules)

plot_gene_heatmap(norm_counts_D9, tfs[1:9])

plot_gene_heatmap(norm_counts_D9, tfs[10:18])

plot_gene_heatmap(norm_counts_D9, chemokines)

plot_gene_heatmap(norm_counts_D30, inh_receptors, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))

plot_gene_heatmap(norm_counts_D30, surface_receptors, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))

plot_gene_heatmap(norm_counts_D30, eff_molecules, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))

plot_gene_heatmap(norm_counts_D30, tfs[1:9], c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))

plot_gene_heatmap(norm_counts_D30, tfs[10:18], c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))

plot_gene_heatmap(norm_counts_D30, chemokines, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))

DGE bw genotypes

select_genes <- c("Btg2", "Ctse", "Cxcr4", "Il17ra", "Tnfrsf9")
dge_D9_progexh %>% my_plot_volcano(label_genes = select_genes) + 
  labs(title = "D9, prog exh")

dge_D9_transexh %>% my_plot_volcano(label_genes = select_genes) + 
  labs(title = "D9, trans exh")

dge_D9_termexh %>% my_plot_volcano(label_genes = select_genes) + 
  labs(title = "D9, term exh")

From single cell too:

select_genes <- c("Btg2", "Ctse", "Cxcr4", "Il17ra", "Tnfrsf9")
sc_dge %>%
  drop_na(p_val_adj) %>% 
  filter(contrast == "D_W") %>% 
  filter(Clusters %in% c("terminal", "progenitor", "transitory")) %>% 
  ggplot() +
  aes(avg_logFC,-log10(p_val_adj)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_hline(
    yintercept = -log10(0.05),
    linetype = "dotted",
    color = "red"
  ) +
  geom_point(aes(color = p_val_adj < 0.05),
             data = ~ ..1 %>% filter(gene %in% c(select_genes)),
             size = 2) +
  geom_text_repel(
    aes(label = gene, color = p_val_adj < 0.05),
    data = ~ ..1 %>% filter(gene %in% c(select_genes)),
    max.overlaps = Inf
  ) +
  scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
  theme(legend.position = "none") + 
  facet_grid(~Clusters) + 
  theme(strip.background = element_blank()) + 
  labs(title = "DGE from scRNAseq")

select_genes <- c("Sell", "Tcf7")
dge_D9_progexh %>% my_plot_volcano(label_genes = select_genes) + 
  labs(title = "D9, prog exh")

dge_D9_transexh %>% my_plot_volcano(label_genes = select_genes) + 
  labs(title = "D9, trans exh")

dge_D9_termexh %>% my_plot_volcano(label_genes = select_genes) + 
  labs(title = "D9, term exh")

From single cell too:

select_genes <- c("Sell", "Tcf7")
sc_dge %>%
  drop_na(p_val_adj) %>% 
  filter(contrast == "D_W") %>% 
  filter(Clusters %in% c("terminal", "progenitor", "transitory")) %>% 
  ggplot() +
  aes(avg_logFC,-log10(p_val_adj)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_hline(
    yintercept = -log10(0.05),
    linetype = "dotted",
    color = "red"
  ) +
  geom_point(aes(color = p_val_adj < 0.05),
             data = ~ ..1 %>% filter(gene %in% c(select_genes)),
             size = 2) +
  geom_text_repel(
    aes(label = gene, color = p_val_adj < 0.05),
    data = ~ ..1 %>% filter(gene %in% c(select_genes)),
    max.overlaps = Inf
  ) +
  scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
  theme(legend.position = "none") + 
  facet_grid(~Clusters) + 
  theme(strip.background = element_blank()) + 
  labs(title = "DGE from scRNAseq")

GSEA

# gsea_res_bp %>%
#   dplyr::select(gs_name, str_subset(colnames(.), "prog")) %>%
#   filter(q.val_progexh < 0.05) %>%
#   arrange(desc(sscore_progexh)) %>%
#   View()
# gsea_res_bp %>%
#   dplyr::select(gs_name, str_subset(colnames(.), "trans")) %>%
#   filter(q.val_transexh < 0.05) %>%
#   arrange(sscore_transexh) %>%
#   View()
dge_list <- list(
  progenitor = dge_D9_progexh %>% mutate(cluster_id = "progexh"), 
  transitory = dge_D9_transexh %>% mutate(cluster_id = "transexh"), 
  terminal = dge_D9_termexh %>% mutate(cluster_id = "termexh")
)
dge_df <- dge_list %>%
  purrr::reduce(bind_rows) %>%
  mutate(signed_neglog10_p = ifelse(log2FoldChange > 0, -log10(pvalue), log10(pvalue))) %>%
  drop_na(log2FoldChange, pvalue)
make_rank_list <- function(.data,
                           gene = gene,
                           rank_by = log2FoldChange) {
  rank_by <- enquo(rank_by)
  gene <- enquo(gene)
  data_fc <- .data %>%
    arrange(desc(!!rank_by))
  data_fc %>%
    .[[rlang::as_name(rank_by)]] %>%
    set_names(data_fc[[rlang::as_name(gene)]])
}
rank_lists <- dge_df %>%
  named_group_split(cluster_id) %>%
  map(~make_rank_list(..1, rank_by = signed_neglog10_p, gene = gene))
my_plot_gsea_curve <- function(gsea_res, rank_list, gs_subset_list, gs_name_str, exh_cluster) {
  my_gs <- gsea_res %>%
    filter(gs_name == gs_name_str)
  Rsc::get_leading_edge_df(rank_vec = rank_list[[exh_cluster]],
                           gs_vec = gs_subset_list[[gs_name_str]]) %>%
    Rsc::plot_gsea_curve(title = paste0(gs_name_str, ", ", exh_cluster)) +
    theme(legend.position = "none") +
    annotate(
      "text",
      x = Inf,
      y = Inf,
      hjust = 1,
      vjust = 1,
      label = paste0(
        "qval=",
        format(my_gs[[paste0("q.val_", exh_cluster)]], digits = 2),
        "\n",
        "sscore=",
        format(my_gs[[paste0("sscore_", exh_cluster)]], digits = 4)
      )
    )
}
my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "REACTOME_CELL_CYCLE", "progexh") + 
    ylim(-0.3, 0.3)

my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "REACTOME_CELL_CYCLE", "transexh") + 
    ylim(-0.3, 0.3)

my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "REACTOME_CELL_CYCLE", "termexh") + 
    ylim(-0.3, 0.3)

my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "progexh") + 
    ylim(-0.3, 0.3)

my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "transexh") + 
    ylim(-0.3, 0.3)

my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "termexh") + 
    ylim(-0.3, 0.3)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_DN", "progexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_UP", "progexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_DN", "transexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_UP", "transexh")  + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_DN", "termexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_UP", "termexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_DN", "progexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_UP", "progexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_DN", "transexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_UP", "transexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_DN", "termexh") + 
  ylim(-0.23, 0.23)

my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_UP", "termexh") + 
  ylim(-0.23, 0.23)

gene_sets <- str_subset(names(gs_subset_bp_list), "RUNX")
gsea_curve_list <- list()
for (gs in gene_sets) {
  gsea_curve_list[[paste0("progexh / ", gs)]] <-
    my_plot_gsea_curve(gsea_D9_bp, rank_lists,
                       gs_subset_bp_list, gs, "progexh")
  gsea_curve_list[[paste0("termexh / ", gs)]] <-
    my_plot_gsea_curve(gsea_D9_bp, rank_lists,
                       gs_subset_bp_list, gs, "termexh")
  
}
gsea_curve_list[c("progexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY", 
                  "termexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY", 
                  "progexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY", 
                  "termexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY", 
                  "progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2",
                  "termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2",
                  "progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1",
                  "termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1")]
## $`progexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY`

## 
## $`termexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY`

## 
## $`progexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY`

## 
## $`termexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY`

## 
## $`progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2`

## 
## $`termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2`

## 
## $`progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1`

## 
## $`termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1`

Rsc::get_leading_edge_df(rank_vec = rank_lists[["progexh"]],
                           gs_vec = gs_subset_bp_list[["REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1"]]) %>% 
  # filter(is_leading_edge) %>% 
  filter(in_gs) %>% 
  write_tsv(here("bulkRNAseq", "results", "leadingedge_progexh_D9_REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1.txt"))
Rsc::get_leading_edge_df(rank_vec = rank_lists[["progexh"]],
                           gs_vec = gs_subset_bp_list[["REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY"]]) %>% 
  # filter(is_leading_edge) %>% 
  filter(in_gs) %>% 
  write_tsv(here("bulkRNAseq", "results", "leadingedge_progexh_D9_REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY.txt"))
Rsc::get_leading_edge_df(rank_vec = rank_lists[["progexh"]],
                           gs_vec = gs_subset_bp_list[["REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY"]]) %>% 
  # filter(is_leading_edge) %>% 
  filter(in_gs) %>% 
  write_tsv(here("bulkRNAseq", "results", "leadingedge_progexh_D9_REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY.txt"))
Rsc::get_leading_edge_df(rank_vec = rank_lists[["termexh"]],
                           gs_vec = gs_subset_tcell_list[["GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_UP"]]) %>% 
  filter(is_leading_edge) %>% 
  filter(in_gs)
## # A tibble: 71 × 7
##    name    value in_gs  rank running_es direction is_leading_edge
##    <chr>   <dbl> <lgl> <int>      <dbl> <chr>     <lgl>          
##  1 Dusp2  -0.501 TRUE  18517     -0.227 dn        TRUE           
##  2 H2ax   -0.512 TRUE  18586     -0.226 dn        TRUE           
##  3 Anln   -0.518 TRUE  18632     -0.223 dn        TRUE           
##  4 Nusap1 -0.542 TRUE  18757     -0.223 dn        TRUE           
##  5 Tmem97 -0.557 TRUE  18824     -0.222 dn        TRUE           
##  6 Ak3    -0.558 TRUE  18830     -0.217 dn        TRUE           
##  7 Dtl    -0.560 TRUE  18839     -0.212 dn        TRUE           
##  8 Tuba1a -0.562 TRUE  18846     -0.207 dn        TRUE           
##  9 Plk4   -0.566 TRUE  18870     -0.203 dn        TRUE           
## 10 Nudt4  -0.567 TRUE  18878     -0.199 dn        TRUE           
## # … with 61 more rows
gs <- "GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_UP"
my_plot_gsea_curve(gsea_D9_tcell, rank_lists,
                   gs_subset_tcell_list, gs, "progexh")

my_plot_gsea_curve(gsea_D9_tcell, rank_lists,
                   gs_subset_tcell_list, gs, "termexh")